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We investigate the rate of thermalization of local operators in the one-dimensional anisotropic 
antiferromagnetic Heisenberg model with next-nearest neighbor interactions that break integrability. 
This is done by calculating the scaling of the difference of the diagonal and canonical thermal 
ensemble values as function of system size, and by directly calculating the time evolution of the 
expectation values of the operators with the Chebyshev polynomial expansion. Spatial and spin 
symmetry is exploited and the Hamiltonian is divided in subsectors according to their symmetry. 
The rate of thermalization depends on the proximity to the integrable limit. When integrability is 
weakly broken thermalization is slow, and becomes faster the stronger the next-nearest neighbor 
interaction is. Three different regimes for the rate of thermalization with respect to the strength of 
the integrability breaking parameter are identified. These are shown to be directly connected with 
the relative strength of the low and higher energy difference off-diagonal operator matrix elements 
in the symmetry eigenbasis of the Hamiltonian. Close to the integrable limit the off-diagonal matrix 
elements peak at higher energies and high frequency fluctuations are important and slow down 
thermalization. Away from the integrable limit a strong low energy peak gradually develops that 
takes over the higher frequency fluctuations and leads to quicker thermalization. 
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I. INTRODUCTION 

The non-equilibrium dynamics of isolated quantum 
systems has recently been an area of extensive investi¬ 
gation m ■ a part of an isolated system feels the rest 
of the system acting as a bath, and is in general in a 
mixed state. Its non-unitary time evolution can lead in 
the long-time limit to a thermalized state where fluctu¬ 
ations of operator expectation values are weak, and this 
state can be described by an appropriate statistical en¬ 
semble depending only on a small number of parameters. 

The conserved quantities of the Hamiltonian are very 
important for thermalization. Classical systems where 
only the energy is conserved are ergodic and thermal- 
ize, and the long time average of local observables can 
be replaced by a time-independent canonical thermal en¬ 
semble average. This describes expectation values in the 
long-time limit and is equivalent to the microcanonical 
ensemble in the thermodynamic limit (TDL). In contrast, 
integrable classical systems are non-ergodic. The number 
of conservation laws scales linearly with system size and 
they confine the motion in phase space to invariant tori 
satisfying the conservation laws. These conserved quan¬ 
tities have to be taken into account in the generalized 
Gibbs ensemble that describes expectation values in the 
long-time limit. When a weak integrabilty breaking term 
is applied the Kolmogorov-Arnold-Moser theorem states 
that some of the invariant tori are deformed and survive 



A formal prescription for the onset of thermalization 
when going away from the integrable limit lacks for quan¬ 
tum systems. The prevailing paradigm is that similarly 
to the classical case non-integrable systems thermalize 
and their long-time limit is described by a canonical 


Gibbs ensemble, while the long-time limit of integrable 
systems is described by a generalized Gibbs ensemble tak¬ 
ing into account all conservation laws 0 03 m El- The 
rate of thermalization is expected to increase with the 
distance from the integrable limit. Questions with re¬ 
spect to thermalization are more challenging to answer 
in quantum mechanics from the analytic as well as the 
numerical point of view, as with present day computa¬ 
tional means it is not possible to calculate the evolu¬ 
tion of observables for very long times. Quite recently 
it has been shown that for weak integrability breaking 
local operator expectation values go through a prether- 
malized regime where they first relax to a non-thermal 
quasi-stationary value before reaching the canonical ther¬ 
mal ensemble value mm- In general, non-equilibrium 
studies of strongly correlated models are confined to a 
relatively small region of parameter space, and very de¬ 
tailed results on the rate of thermalization as function of 
integrability breaking parameters are not available. 

The long-time limit of an operator expectation value 
is given by the operator’s diagonal matrix elements in 
the energy eigenbasis, which define the diagonal ensemble 
mini. Thermalization requires the difference between 
the diagonal and the canonical Gibbs ensemble to vanish 
in the TDL. The strength of the fluctuations around the 
long-time limit is controlled by the off-diagonal terms 
of the operator matrix, which determine if fluctuations 
are small enough in order to achieve thermalization, and 
also which time scales are important for relaxation. Re¬ 
sults on this topic have been presented in the literature 
pn US], and quite recently the statistical properties of 
the off-diagonal matrix elements have been studied for 
two models with respect to their proximity to integrabil¬ 
ity m- However, these authors intentionally break spa- 
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tial symmetry and consequently their results are neither 
representative of typical connectivities of strongly cor¬ 
related models nor of individual symmetry sectors that 
oftentimes solely include the time-evolved wavefunction. 

In this paper we investigate the rate of thermalization 
when integrability breaking is controlled by a parameter. 
We show how thermalization of local operators, which 
relates to the time dependence of the fluctuations of the 
operator expectation values, is directly connected to the 
strength of the operator off-diagonal matrix elements. 
We do so within the framework of the anisotropic antifer¬ 
romagnetic Heisenberg model (AAHM) on chains with N 
spins of magnitude s = 1/2 and periodic boundary con¬ 
ditions. The model is taken to have nearest neighbor 
interactions, and weaker next-nearest neighbor interac¬ 
tions that break integrability in a controlled manner by 
tuning their strength. We consider the typical in non¬ 
equilibrium dynamics case of an interacting quench and 
take the spatial and spin symmetries of the Hamiltonian 
into account. The initial state is the ground state of the 
AAHM for specific anisotropy, chosen so that the post- 
quench effective temperature is approximately the same 
for all values of the integrability breaking parameter of 
the model. The Hamiltonian is diagonalized according to 
its full symmetry group before and after the quench [18}- 
[20] , and the difference of the diagonal and thermal en¬ 
semble values for local operators is calculated as function 
of system size to establish if it tends to zero at the TDL, 
as required by thermalization. In addition, the time evo¬ 
lution of the operators is calculated with the Chebyslrev 
polynomial expansion mm for N significantly larger 
than the system maximum that can be exactly diago¬ 
nalized and for long enough times that show how close 
the system is to thermalization. The results show how 
the rate of thermalization increases with the integrability 
breaking parameter. We also calculate the expectation 
values of the off-diagonal elements of the operators in 
question in the symmetry eigenbasis of the Hamiltonian, 
and more specifically within the symmetry subsector that 
exclusively includes the wavefunction after the quench. 
Their distribution as function of the energy difference of 
the eigenstates is shown to be in direct relation with the 
rate of thermalization. The observables considered are 
the spin correlation functions sfsj + sfsj and sfsj that 
are a maximum of two sites apart \j — i\ = 1,2 along the 
closed chain. 

We find that for weak integrability breaking the differ¬ 
ence between the diagonal and canonical thermal ensem¬ 
ble values decreases sublinearly with indicating slow 
thermalization. Here it is assumed that the equilibrium 
value of an operator can be found by first calculating the 
long time average and then taking the TDL [23], which 
is expected to be adequate for large interacting systems. 
The corresponding time-dependent expectation values of 
the operators in question do not differ significantly from 
the thermal expectation values, but these differences do 
not decrease significantly for the highest available times 
which further corroborates a relative slow convergence 


to the thermal ensemble for weak integrability break¬ 
ing. When the integrability breaking term gets stronger, 
roughly one tenth of the nearest neighbor interaction, 
thermalization becomes faster with the scaling of the dif¬ 
ference between the two ensembles very close to linear. 
Furthermore, the time evolution of the operator expec¬ 
tation values approaches the thermal ensemble with a 
faster pace. For even stronger integrability breaking the 
scaling of the difference of the two ensembles becomes su- 
perlinear and the time dependent operators approach the 
thermal ensemble values very closely. We conclude that 
there are three different regimes with respect to the rate 
of thermalization when integrability is broken, which are 
in direct correlation with the strength of the integrabil¬ 
ity breaking parameter. This conclusion is supported by 
two different calculations that are in very close agreement 
with each other. 

The rate of thermalization is determined by the 
strength of the operator off-diagonal elements in the en¬ 
ergy eigenbasis, which determine the magnitude of tem¬ 
poral fluctuations. When operator matrix elements be¬ 
tween states that differ significantly in energy strongly 
prevail over their counterparts with small energy differ¬ 
ences, high frequency fluctuations are strong and ther¬ 
malization is slow. The reason is that it takes a consid¬ 
erable time for these fast fluctuations to cancel each other 
out and be overtaken by the weak low frequency fluctua¬ 
tions that correspond to matrix elements between eigen¬ 
states with small energy difference. In contrast, if matrix 
elements between states not differing much in energy are 
also of considerable strength high frequency fluctuations 
diminish faster, leading to quicker thermalization. To 
explain the thermalization scenario emerging from the 
comparison between the difference of the diagonal and 
thermal ensembles and the time evolution of operator 
expectation values we calculate the distribution of the 
corresponding operator off-diagonal matrix elements as 
function of energy difference in the symmetrized energy 
eigenbasis. In the integrable limit the Hamiltonian sub¬ 
sector where time evolution takes place is characterized 
by the required spatial, total spin along the z axis S z and 
spin inversion symmetries (SU(2) symmetry is in general 
broken for the AAHM) of the Hamiltonian, however there 
are more conserved quantities 1241 which are not taken 
into account and consequently do not characterize the 
states in the subsector of the time evolved wavefunction. 
A Hamiltonian subsector where all symmetries have been 
taken into account is characterized by level repulsion and 
has an energy difference between any pair of states which 
is in general not close to zero, provided there are no 
Hamiltonian parameters much smaller than the others. 
On the other hand, in the integrable limit where many of 
the conservation laws are not taken into account in the 
block diagonalization of the Hamiltonian it is expected 
that there will be levels very close in energy, as they 
would belong to different subsectors of the Hamiltonian if 
the extra conservation laws were taken into account. Be¬ 
cause of that the Hamiltonian matrix elements between 
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such levels should be close to zero. The nearest neigh¬ 
bor correlation functions sfsf, 1 + s^sf +1 and s?sf +1 are 
translationally invariant in singly degenerate subsectors 
and qualitatively similar to the nearest neighbor energy, 
consequently their off-diagonal matrix elements between 
eigenstates very close in energy should also be close to 
zero in the integrable limit. Hence the off-diagonal ele¬ 
ments of these operators have a peak for an energy dif¬ 
ference away from zero. 

Once the integrability breaking term is switched on 
the additional conservation laws not taken into account 
in the integrable case are not fulfilled anymore, conse¬ 
quently states with small energy differences will start 
mixing more strongly and the correlation functions will 
develop significant off-diagonal elements between such 
states. The strength of these off-diagonal elements is 
regulated by the magnitude of the integrability break¬ 
ing term. According to low order perturbation theory 
the mixing is in fact favored for levels which are close in 
energy in the integrable limit. This results in a peak for 
small energy differences that now coexists with the higher 
energy peak that determines the higher frequency fluc¬ 
tuations. When integrability is weakly broken the low 
energy peak is relatively weak in comparison with the 
higher energy peak. This shows that a long time is re¬ 
quired for the high frequency fluctuations to cancel each 
other out and thermalization to settle in, resulting in a 
long transient regime before thermalization. For stronger 
next-nearest neighbor interaction the low energy peak 
becomes of the same order of magnitude with the higher 
energy peak, thus the high frequency fluctuations become 
unimportant faster and thermalization is quicker. This 
provides the explanation to the results calculated with 
exact diagonalization and the Chebyshev polynomial ex¬ 
pansion. All in all, the non-equilibrium properties of the 
AAHM with a next-nearest neighbor integrability break¬ 
ing term are determined by the strength of this term. 

The plan of this paper is as follows: In Sec. |TT| the 
model and the methods to be used are introduced. In 
Sec. |III| general considerations about the time evolution 
of an operator expectation value after a quench are pre¬ 
sented. Sec. |IV| examines the scaling of the difference of 
the diagonal and canonical ensemble for an interacting 
quench and also includes time-dependent operator ex¬ 
pectation values. Sec. [V] correlates thermalization with 
the strength of operator off-diagonal terms, and finally 
Sec. \VJ\ presents the conclusions. 


II. MODEL AND METHOD 

The Hamiltonian of the AAHM model on a chain with 
N spins is 

N 

H = ^ Sj * S?+l "F J ‘ Si-1-2) ( 1 ) 

i=1 


The spin magnitude is taken to be s = 1/2 and periodic 
boundary conditions are assumed so that Sat+i = Si, 
sn +2 = S 2 . N is taken to be even. J controls the strength 
of nearest neighbor interactions, and J' the strength of 
next-nearest neighbor interactions that break integrabil¬ 
ity. Both of them are taken positive and thus support 
antiferromagnetic correlations, with J' introducing frus¬ 
tration. The interaction along the spin z axis is scaled 
by A, Si ■ Sj = sfsj + sfsj + A s z s z . J = 1 from now on, 
defining the unit of energy, and A' = A. 

Calculation of the diagonal and thermal ensemble val¬ 
ues requires full diagonalization of Hamiltonian Q. To 
this end, the spatial Dn and the spin inversion symme¬ 
try (when S z = 0) are exploited, which is particularly 
advantageous [T3H20] . The time evolution of various op¬ 
erators can be calculated directly with exact diagonaliza¬ 
tion, but the more economical method of the Chebyshev 
polynomial expansion increases the maximum N that can 
be considered especially since the Hamiltonian 

block-diagonalizes according to the irreducible represen¬ 
tations of the total symmetry group and the wavefunction 
is time-evolved separately within each representation. If 
the initial wavefunction belongs solely to a specific sym¬ 
metry subsector, then it time evolves only in this sub¬ 
sector and also the diagonal ensemble contains operator 
expectation values of this subsector only. This is the case 
here as the initial state in an interacting quench is the 
ground state of the initial Hamiltonian and therefore be¬ 
longs to a specific irreducible representation of the full 
symmetry group of Hamiltonian 0 - The initial state 
is chosen to be the ground state of Hamiltonian |T]) for 
Ao = 15, chosen so that the post-quench effective tem¬ 
perature is nearly the same for all values of J' , with the 
post-quench A = 1.1. 

Hamiltonian 0 is equivalent to the one of the saw¬ 
tooth chain with variable interactions fToWIEl . If J' = 0 
or —> oo it corresponds to the unfrustrated AAHM 
and is integrable via the Bethe Ansatz [29]. For A = 1 it 
is also integrable for the Majumdar-Glrosh point — \ 
[301 [3Tj . For A = 1 and small J j there is quasi-long range 
order in the ground state which gives way to a gapped 
phase characterized by a dimerization order parameter 
at j slightly less than | [32]. Here we focus on the re¬ 
gion of relatively small iy < | to investigate the effect 
of integrability breaking in thermalization for operators 
of the form sfsj + sfs^ and sf s| with \j — i\ = 1, 2. We 
pick A = 1.1 for the post-quench Hamiltonian to avoid 
the mixing of the different S sectors of the spin isotropic 
SU(2) case when one works in the S z basis. 


III. TIME EVOLUTION AFTER A QUENCH 

The initial wavefunction |\E'(t = 0)) = |'F(0)) is chosen 
to be the ground state before the quench. Its post-quench 
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time evolution is given by 

|*(t)> = e~ iHt \^(0)) =Y J C n e~ iE - t \^ n ) ( 2 ) 

n 

with E n and |4/ n ) the eigenvalues and eigenvectors of 
Hamiltonian (|T]) after the quench. The coefficients C n = 
(4'„|\E'(0)) are independent of time. The expectation 
value of operator O at time t is 

< 0(t) > = £|C'„| 2 <*„|0|* J1 ) + 

n 

Y |0|$ m ) (3) 

m^n 

In the time average of Eq. © the second term 
averages out to zero and the average is O = 
J2n |Cn| 2 (^ / n| 0 |’I'n), defining the diagonal ensemble 
Odiag = ^„(^ , n|0|^'n)|^ , n)(^ , n| [H] In the interacting 
quench considered here time evolution takes place in non¬ 
degenerate subsectors, therefore degeneracies play no role 
in the second term of Eq. 0. For thermalization to 
occur the diagonal ensemble prediction must equal the 
value of the canonical thermal ensemble in the TDL. Eq. 
([3]) shows that the strength of the fluctuations around 
the diagonal ensemble value is determined by the off- 
diagonal terms (4'„|0| v h m ) = O nm , m ^ n. Another 
necessary condition for thermalization is that the fluctu¬ 
ations with respect to 0 must be small in the TDL. The 
prevailing paradigm for the thermalization mechanism is 
the Eigenstate Thermalization Hypothesis [3], but more 
general ideas have also been proposed [23] ■ To calculate 
the canonical thermal ensemble p t h that describes the 
equilibrated long time limit the energy (\E r (0)|.ff|\I r (0)) is 
set equal to the thermal average Tr(Hp t h), and the equa¬ 
tion is numerically solved for the temperature. 

IV. SCALING AND TIME DEPENDENCE OF 
LOCAL CORRELATIONS 

To investigate how the rate of thermalization depends 
on the integrability breaking parameter J' of Hamilto¬ 
nian 0 the difference of the diagonal and canonical 
thermal ensemble values for a local operator O is cal¬ 
culated as function of the size of the chain N. More 
specifically, the relative difference of the two ensembles 
is defined as the difference in their values divided by 
the value of the canonical thermal ensemble, so that 
A O re i = - The scaling of A(sfsf +1 + sfsf +1 ) re ; 

with A is shown in Fig. 0a) as function of J' for N < 24. 
A(sf sf +1 +sfsi +1 ) re i is getting smaller in magnitude with 
increasing N and appears to be directed towards the ori¬ 
gin in the TDL where —> 0, which agrees with the 
expectation that the Hamiltonian of Eq. 0 thermalizes 
as long as J' 7 ^ 0. If a power law dependence ajAf) b is 
fitted through the data for each J' value, Fig. [2] shows 
the dependence of the exponent b on J'. For very small 


J' it is b < 1 indicating that convergence to the thermal 
ensemble is slow for weak integrability breaking. As J' 
increases b becomes approximately 1 around J' ~ 0.1, 
showing that convergence becomes linear and thermal¬ 
ization faster. For even stronger J' it is b > 1 and ther¬ 
malization becomes even faster. 

These results can be corroborated by direct calcula¬ 
tion of the time evolution of < sfsf +1 + s^s v i+1 > with 
the Chebyshev polynomial expansion mm which can 
treat larger chains in comparison with exact diagonaliza- 
tion. This data is shown in Fig. [3] and the difference 
between the time evolved < sfsf +1 + sfs ) / +1 > and the 
canonical thermal ensemble value scaled with the canon¬ 
ical thermal ensemble value is shown in Fig. [4] The 
finite size of the system eventually generates revivals of 
the operator expectation value in later times, and the 
closed chain with N = 30 used here generates the infi¬ 
nite chain results for times at least up to ® ; 12] . This has 
also been vindicated by comparison of the results as func¬ 
tion of V, which agree to progressively longer times with 
increasing N. The slow thermalization close to the inte- 
grable limit has brought forward the concept of prether- 
malization, where local operator expectation values go 
through a prethermalized regime where they first relax 
to a non-thermal quasi-stationary value before reaching 
the canonical thermal ensemble value mm- This does 
not imply absence of thermalization, but rather that fluc¬ 
tuations are more slowly varying in time in comparison 
with stronger integrability breaking. Here it is also ob¬ 
served that for weak J' < 0.05 (Figs[3]and0 there is slow 
convergence of the time expectation value to the thermal 
ensemble, which agrees with the results for the scaling of 
A(s? sf +1 +s\ S V i+1 )rel with that gave an exponent b < 1 
in the power law fit. For stronger J' = 0.1 where the fi¬ 
nite size scaling of A(sfsL 1 + s^s^ +1 ) re i points to faster 
thermalization Figs [3] and [4] also show that, with the time 
dependent difference tending to zero more quickly and 
the smaller time fluctuations becoming weaker. For the 
highest values of J' Figs[3]and[4]show an even higher ther¬ 
malization rate, in agreement with the scaling data. The 
difference between the time evolved < sfsf +1 + s\s v i+1 > 
and the thermal ensemble is becoming smaller with an 
even faster rate. 

A similar procedure for the nearest neighbor correla¬ 
tion function S i S i +1 generates the scaling data of Fig. 
0b). Again A(s|sf +1 ) re i decreases with N and ap¬ 
pears to be directed towards the origin in the TDL 
where —> 0. A power law fit a(jj) b shows that the 
rate of thermalization changes as function of integrabil¬ 
ity breaking in exactly the same way as with correlation 
s i s f+i + s f s f+i> the exponent b having the same 

dependence on J' and being close in value to the one of 
A(sfsL i + sfs® , + i)r-ei (Fig.©. The time-dependent data 
of Fig. ©also show that for the maximum available time 
the difference between the expectation value < s£ sf +1 > 
and the canonical thermal ensemble value decreases with 
J'. 

Considering next-nearest neighbor correlations, the de- 
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pendence of A(sfsf +2 + sfsf +2 ) re z and A(sfsf +2 )„.j on 
T is shown in Fig. [6j They decrease in magnitude as 
N increases, and the exponent b of a fit to a function 
of the form a(-^) b is plotted in Fig. rd as function of 

. Again three different regimes witn respect to the 
rate of thermalization as function of N are seen, with 
the exponents for the two operators not differing much. 
The time-dependent data of Figs [7] and [8] confirm the 
similarity with the nearest neighbor operator case, with 
the difference between the time-dependent expectation 
values and the canonical thermal ensemble value becom¬ 
ing smaller when going away from the integrable limit as 
thermalization becomes faster. 

Comparing the discrepancy of the time expectation 
values from the canonical thermal emsemble value, it 
is seen from Figs [3] and [7] and Figs [5] and [8] that the 
nearest neighbor correlations approach their thermalized 
limit faster than the next-nearest neighbor correlations 
for the longest times available. This is because nearest 
neighbor correlations are more local, therefore the rest of 
the system acts as a bigger and more effective bath, at 
least for the sizes available for the calculations. 

Comparing the value of the fitting exponent b between 
the two cases in Fig. [2j for weak J' the exponent is 
smaller for next-nearest neighbor correlations. This has 
to do again with the more local nature of nearest neigh¬ 
bor correlations. Since convergence to the TDL limit is 
slow, it will take a larger N before the effective bath size 
becomes equivalent for correlations betweens spins one 
and two sites way. With increasing J' the exponent b be¬ 
comes higher for next-nearest neighbor correlations. For 
higher J' the interactions are becoming less local and J' 
directly connects spins two sites apart in a more efficient 
way. Thus progressively a smaller system size is required 
for the effective bath size to be the same for two and three 
sites, consequently b is higher for next-nearest neighbors 
as they approach this system size fast on the way to the 
TDL. 


V. OFF-DIAGONAL TERMS AND 
THERMALIZATION 


As discussed in 


Sec. Ill Eq. ^ demonstrates that the 


off-diagonal terms O nm determine the rate of equilibra¬ 
tion and the strength of fluctuations around the infinite 
time limit. Quantum chaos theory predicts that the O nm 
are characterized by a smooth distribution [13], [14] . More 
precisely, the fluctuations in Eq. are determined by 
the strength of the O nm and the phases E n — E m = 
A E nm . High frequency fluctuations involve eigenstates 
with large A E nm . At longer times these terms average 
out to zero and the long time limit is determined by over¬ 
laps between eigenstates close in energy. The relative 
strength of O nm 's with small and larger energy differ¬ 
ences is therefore crucial for the rate of equilibration. 
The distribution of the (sfsf+i + s i s i+i)nm with A E nm 
depends on the subsector of the Hamiltonian that time 


evolution takes place and if all symmetries have been 
taken into account in the time evolution, which is de¬ 
termined by the integrability or not of Hamiltonian 
and also on the specific values of the Hamiltonian pa¬ 
rameters, which determine the strength of the overlaps 
between different eigenstates. 

The Hamiltonian subsector where time evolution takes 
place is determined by the ground state before the 
quench. This subsector belongs to a specific one- 
dimensional irreducible representation of the total sym¬ 
metry group of Hamiltonian |l]). However, in the in¬ 
tegrable limit there is a total number of N conserved 
quantities which are not fully taken into account. Conse¬ 
quently the Hamiltonian subsector that determines time 
evolution has states that would belong to different sub¬ 
sectors had all the conservation laws been taken into 
account. The nearest neighbor correlation functions 
sfsf +1 + s v i s v i+l and s* sf +1 are translationally invariant 
in singly degenerate subsectors and qualitatively simi¬ 
lar to the nearest neighbor energy (it is noted that their 
sum does not correspond to the energy when A / 1). A 
Hamiltonian subsector where all symmetries have been 
taken into account is characterized by level repulsion and 
has an energy difference between any pair of states which 
is in general not close to zero, provided there are no small 
parameters in the Hamiltonian. As in the integrable limit 
many of the conservation laws are not taken into account 
in the block diagonalization of the Hamiltonian it is ex¬ 
pected that there will be levels very close in energy, which 
would belong to different subsectors of the Hamiltonian if 
the extra conservation laws were to be taken into account. 
Because of this reason the Hamiltonian matrix elements 
between such levels should be close to zero, and conse¬ 
quently the same should hold for the off-diagonal matrix 
elements of correlations sfsf +1 + sfsf +1 and sfsf +1 . For 
this reason in the integrable limit the distribution of the 
off-diagonal elements of the spin correlation functions has 
small values close to zero and peaks at an energy differ¬ 
ence away from zero. 

Once integrability is broken with J' ^ 0 the additional 
conservation laws of the integrable case do not apply 
anymore. States with small energy differences start to 
mix more strongly and the correlation functions start to 
develop significant off-diagonal matrix elements between 
such states, especially since according to low order per¬ 
turbation theory the mixing is stronger for levels having 
close energies in the integrable limit. The strength of 
these off-diagonal elements is regulated by the magni¬ 
tude of the integrability breaking term. A plot of the 
(sfsf j_l + sfs* / +1 ) nm as function of A E nm is shown in 
Fig. |9]for N = 18 for J' = 0.01, 0.1 and 0.2, along with 
a coarse grained average of the matrix elements as func¬ 
tion of A E nm in Fig. [9] (d). For weak J' = 0.01 the 
corresponding distribution of Fig. |9^a) and the coarse 
grained average in Fig. [9)(d) peak away from the origin 
as function of A E nm , with a much weaker low energy 
peak induced by the integrability breaking term. Ac¬ 
cording to Eq. ([3]) since the high-frequency fluctuations 
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have a significantly higher strength it takes a consider¬ 
able time for them to sufficiently diminish due to the 
mixing of their different phases and only then the signif¬ 
icantly weaker low frequency modes dominate time evo¬ 
lution. For the short times available to numerical cal¬ 
culations it is not possible to see later times where the 
long time limit has settled in, but only a transient regime 
where the system has not yet fully thermalized and the 
expectation value has not fully converged to its infinite 
time limit (Figs [3] and [4]). A bigger J' = 0.1 generates 
stronger (sfs? +1 + s* , s® / +1 )„ m ’s for low energy differences 
and for even higher A E nm as evidenced by Figs[9](b) and 
[9jd). The strength of the low energy peak becomes of 
the same order of magnitude with the high energy peak 
and thermalization is now faster. The effects become 
even stronger for J' = 0.2. The low-energy peak acquires 
strength comparable to the higher energy peak and the 
two peaks also begin to merge (Figs |9](c) and (d)). This 
results in the qualitative change in the time evolution of 
< sfsf +1 + s‘ / sf +1 > and its faster thermalization. The 
above provide the explanation for the scaling of the cor¬ 
relations functions and the time-dependent data. It must 
be noted that the above considerations do not take into 
account the energy distribution of the coefficients of the 
initial state C n (Eq. ([3])), however their distribution has 
bigger weight for smaller A E nm as they approach a sharp 
peak going to the TDL, therefore at most renormalize the 
low energy peak of the (sfsf +1 + sfsf +1 )„ m distribution 
and strengthen it with respect to the higher energy peak. 

VI. CONCLUSIONS 

The approach to thermalization as function of integra- 
bility breaking was investigated for the one-dimensional 


AAHM with nearest and next-nearest neighbor interac¬ 
tions, with spatial and spin symmetries taken into ac¬ 
count. The scaling of the difference of the diagonal and 
canonical thermal ensembles was examined for local op¬ 
erators, whose time evolution was also calculated for even 
larger systems. Both pointed to a rate of thermalization 
depending on the strength of the next-nearest neighbor 
interactions J'. Three different regimes were identified, 
the first for weak J' where thermalization is slow and the 
dependence of the difference of the two ensembles on 
is sublinear, then a second for stronger J' where ther¬ 
malization is faster and the difference is approximately 
linear in j ?, and then a third where thermalization is even 
faster for even larger J' with superlinear dependence of 
the difference on A. The existence of three different ther¬ 
malization regimes was explained by the dependence of 
the off-diagonal elements of the operators in question on 
the energy difference of the corresponding eigenstates. 
They were shown to develop a low energy peak that de¬ 
termines the long time fluctuations, grows stronger with 
J' and competes with the high energy peak facilitating 
faster thermalization as its strength increases. 


The author is grateful to A. Lazarides for many sug¬ 
gestions and comments, and thanks D. Schuricht for a 
discussion. 
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FIG. 1: Scaling of the difference A O re i = — ^' >th between 
the diagonal and canonical thermal ensemble values divided 
by the canonical thermal ensemble value as function of inverse 
length A f or correlation function (a) O = s^sf +1 + s y i s v i , 1 and 
(b) O = .sf.s.f +1 for A = 1.1 and different values of J’. The 
initial state is the ground state for Ao = 15. The maximum 
chain has N = 24. 
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FIG. 2: Exponent & of a power law A O re i = a(-^) 6 fitted 
through A Orel = as function of J' with A = 1.1 

0 __ \(-')th 

and Im . The initial state is the ground state for 

Ac = 15. 



FIG. 3: Time evolution of < sfsf +1 + s y s y i+1 > for dif¬ 
ferent values of J' with A = 1.1. The dashed lines cor¬ 
respond to the values of the canonical thermal ensemble 
< sfsf+i + s y s v i+1 >th- The initial state is the ground state 
for Aq = 15. 



















FIG. 4: Time evolution of the expectation value of O = 
sfsf +1 + s“sV + j with respect to its thermal value < O >th 
for different values of J' with A = 1.1. The initial state is the 
ground state for Aq = 15. 



FIG. 5: Time evolution of < sfsf +1 > for different values of 
J' with A = 1.1. The dashed lines correspond to the values 
of the canonical thermal ensemble < sfsf +l > t h- The initial 
state is the ground state for Ao = 15. 
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FIG. 6: Scaling of the difference A O re i = °^' >th between 
the diagonal and canonical thermal ensemble values divided 
by the canonical thermal ensemble value as function of inverse 
length -b for correlation function (a) O = sf sf +2 + s^ +2 an( i 

(b) O = sfsf +2 for A = 1.1 and different values of J’ . The 
initial state is the ground state for Ao = 15. The maximum 
chain has N = 24. 



FIG. 7: Time evolution of < sfsf +2 + s^sf +2 > for dif¬ 
ferent values of J' with A = 1.1. The dashed lines cor¬ 
respond to the values of the canonical thermal ensemble 
< sfsf +2 + sfsf + 2 > t h- The initial state is the ground state 
for Ao = 15. 
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FIG. 8: Time evolution of < sfsf +2 > for different values of 
.J' with A = 1.1. The dashed lines correspond to the values 
of the canonical thermal ensemble < sfsf +2 >th- The initial 
state is the ground state for Ao = 15. 



FIG. 9: Distribution of matrix elements (s^sf +1 + sfsL_ 1 ) n m 
as function of the energy difference A E nm = E n — E m for 
(a) J' = 0.01, (b) J' = 0.1, and (c) J' = 0.2. The distri¬ 
bution is symmetric with respect to the A E nm axis. The 
coarse grained average of + sfSi +1 )nm | divided with 

its maximum value is plotted in (d) as function of the reduced 
average energy difference for 100 bins with an equal number 
of points. 






























